function[Results,RobustInter,RobustPars,RobustValues,RobustNE,NPLiterations,Lemondiff] = estimation_func(data,PAR)
% Estimation function
% Input: estimation options (PAR), data, which sample
% Output: structural parameter estimates, NPL diagnostics and convergence, Lemon effect
% ----------------------------

 try
     rng('default') % Fix random seed in estimation for random noise starting values thetas
         if PAR.highvalue == 0 % Price bracket 0-200 AND nov 2022 added reserve price <200 / main sample =======================
                data      = data(data.Sample_main==1,:);
                PAR.th0   = PAR.th0_main;
                % Fees
                cB = 0 ;
                cS = 0.085*1.2;  % For homogenized values: seller commission equal to the [0,200] value bracket, 8.5% + VAT
                cL = 1.75*1.2;
                cR = 0.50*1.2;   % Charging 25 pence for secret reserve and 50 for min bid, most have secret reserve.
                PAR.trimquality=PAR.trimQmain;
                PAR.truncu0hat=PAR.truncu0main;

         else                 % Price bracket 200-1500 AND nov 2022 added reserve price <800 / high value sample =======================
                data      = data(data.Sample_highend==1,:);
                PAR.th0   = PAR.th0_high;
                % Fees
                cB = 0 ;
                cS = 0.075*1.2;  % For homogenized values: seller commission equal to the (200-800 (until 1500)]= 7.5%+VAT
                cL = 1.75*1.2;
                cR = 0.50*1.2;   % Charging 50 pence for secret reserve and 25 for min bid.
                PAR.trimquality=PAR.trimQhigh;
                PAR.truncu0hat=PAR.truncu0high;
        end

        % Prepare for first-stage regression ("homogenization step")
        Tq = (data.numberbidders>1)&((data.highestbid~=data.reserve)|((data.highestbid==data.reserve)&data.sold==0));         % -- Select relevant auctions for estimation of gz
        data.log_highestbid_bot = log(data.highestbid_bot);
            
        % Including N dummies
        data.N0 = data.numberbidders==0;
        data.N1 = data.numberbidders==1;
        data.N2 = data.numberbidders==2;
        data.N3 = data.numberbidders==3;
        data.N4 = data.numberbidders==4;
        data.N5 = data.numberbidders==5;
        data.N6 = data.numberbidders==6;
        data.N7 = data.numberbidders==7;
        data.N8 = data.numberbidders==8;
        data.N9 = data.numberbidders==9;
        data.N10 = data.numberbidders==10;
        data.N11 = data.numberbidders==11;
        data.N12 = data.numberbidders==12;
        data.N13 = data.numberbidders==13;       
          
        % Run OLS regression
        PAR.Firststage_model2=[PAR.Firststage_model,'+N3+N4+N5+N6+N7+N8+N9+N10+N11+N12+N13'];
        LM      = fitlm(data(Tq,:),PAR.Firststage_model2);
        LMcoef  = LM.Coefficients.Estimate;
  
        % Get g(Z) and residual Values
        data.Quality = predict(LM,data)-LMcoef(1);            % Estimated g(Z)
        data.Uhat    = log(data.highestbid_bot)-data.Quality; % Estimated conditional valuation: log(V)=g(Z)+U

        % Drop rows with extreme estimates
        data      = data(quantile(data.Quality,PAR.trimquality)<=data.Quality&quantile(data.Quality,1-PAR.trimquality)>=data.Quality,:);
        
        % Select positive reserve sample
        data_posr = data(data.hasreserve==1,:);

        % Lambdahat(r=0) : ML / simply the first moment if Poisson, otherwise fitting Normal distribution
        lambdastarnor   = mean(data.numberbidders(data.hasreserve==0));
           
        % Bidder's conditional value distribution : MLE 
        thetab = fmincon(@(theta) bidderLL(theta,lambdastarnor,data,PAR), PAR.th0,[],[],[],[],PAR.lb,PAR.ub,[]);

        % Lambdahat(positive reserve (>1)) : MLE
        parameters  = fmincon(@(parameters) bidderLLlambda(thetab,parameters,data_posr,PAR), [lambdastarnor-3 0],[],[],[],[],[0.1 0],[lambdastarnor+5 1],[]);
        lambdastarposr = parameters(1);     %first parameter generalized Poisson
        PAR.sharenobidders = parameters(2); %second parameter generalized Poisson
         
        % Implied conditional seller values in positive reserve auctions
        rhomtilde = log(data_posr.reserve_bot)-data_posr.Quality;
        ggdCDF = arrayfun(@(x) ggd_cdf(x,thetab),rhomtilde);
        ggdPDF = arrayfun(@(x) ggd_pdf(x,thetab),rhomtilde);
        markup = (1-ggdCDF)./ (((1./exp(data_posr.Quality))*(1+cB).*ggdPDF));%returns NAN for rhomtilde outside of support GGD where PDF=0
           
        % lnV0>0 so rpos-markup cannot be negative. Exclude values from analysis (will affect model fit)
        keepr     = ((data_posr.reserve_bot-markup)>0&(imag(markup)==0)&~isnan(markup)&abs(markup)<Inf);
	    data_posr = data_posr(keepr,:);
        deleted_negative_impliedlnV0 = (1-sum(keepr)/length(keepr)); 
        data_posr.log_residualr = log((data_posr.reserve_bot-markup(keepr))*(1-cS));
        data_posr.Quality_seller = data_posr.Quality;
        data_posr.U0hat = data_posr.log_residualr-data_posr.Quality_seller;
	
	    % Get screening value and initial seller entry threshold
        u0star_init    = quantile(data_posr.U0hat,1-PAR.truncu0hat);
        screeningvalue = quantile(data_posr.U0hat,PAR.truncu0hat);

        % Estimate initial seller parameters using likelihood function concentrated at u0star_init
        % sensitive to starting points, especially with Generalized
        % Gaussian distribution. Multistart 1000 startpoints, stop when no error. 
        % Otherwise just an unlucky draw of small sample (esp. bootstrap) :-- quit this run and go to next
	    npmean=mean(data_posr.U0hat);
		npsd=std(data_posr.U0hat);
	    startpts      = zeros(1000,3);
        normmat  = random('normal',npmean,npsd,1000,1000);
	    startpts(:,1) = normmat(:,1)+random('unif',-2,2,1000,1);
        for i = 1:1000
		    startpts(i,1)=max(0,startpts(i,1));
	    end
	    startpts(:,2) = std(normmat,[],2)+random('unif',-2,2,1000,1);
	    for i = 1:1000
		    startpts(i,2)=max(0.05,startpts(i,2));
	    end
        startpts(:,3) = random('unif',-0.5,0.5,1000,1);
        	    
	    count = 1;
        err_count = 1;
        differences=zeros(5,1);
        ub=[8 4 2];%initiate some huge values
        thetas_init = ub;
        while (err_count == 1 ||  sum(thetas_init<ub)<length(thetas_init)) && count < 1001
            try
               [thetas_init,fval] = fmincon(@(theta) sellerLL(theta,PAR,data_posr,u0star_init,screeningvalue), startpts(count,:), [],[],[],[],PAR.lb,PAR.ub,[]);
               err_count = 0;
            catch
            end
            count = count + 1;
        end
        CDFu0star_init = ggd_cdf(u0star_init,thetas_init);
     	
        % Next: Simulate listing-level bidder and seller values given different number of bidders per listing and levels of v0
        % Pre-calculate grid of bidder and seller expected surplus for (n,v0)
        [pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec]  = piBS_MC(cB,cS,thetab,thetas_init,PAR);

        % Estimate opportunity cost of time bidders
        eBo_posr_init = eBofunc(screeningvalue,u0star_init,lambdastarposr,pibmat,pib0vec,pdfmat,v0vec,0,PAR); % second-last input: zeroreserve=0
        eBo_nor  = eBofunc(screeningvalue,u0star_init,lambdastarnor,pibmat,pib0vec,pdfmat,v0vec,1,PAR);  % second-last input: zeroreserve=1       
            
        % Initial lambdastar_posr and eSo
        lambdastarposr_init = lambdastarposr;
        eSo_init      = max([eBo_posr_init,eBo_nor]); 
      
        % NPL iteration or update v0* just once
        tol = PAR.NPLtol;
        maxit = PAR.NPLmaxit;            
        diff = 2;
        iter = 1;
        
        % Initial values
        eSo=eSo_init;
        NPLiterations = zeros(maxit+1,9);
        CDFu0star = CDFu0star_init;
        output = [iter-1,CDFu0star_init,u0star_init,lambdastarposr_init,eSo_init,thetas_init,fval];
        NPLiterations(iter,:)=output;
        thetas=thetas_init;

	    startpts      = zeros(1000,3);
        normmat  = random('normal',mean(thetas(1)),thetas(2),1000,1000);
	    startpts(:,1) = normmat(:,1)+random('unif',-2,2,1000,1);
        for i = 1:1000
                startpts(i,1)=max(0,startpts(i,1));
        end
        startpts(:,2) = std(normmat,[],2)+random('unif',-2,2,1000,1);
        for i = 1:1000
                startpts(i,2)=max(0.05,startpts(i,2));
        end
        startpts(:,3) = random('unif',-0.5,0.5,1000,1);
        startpts(1,:)=thetas;
        
        % Run NPL until convergence
        while diff > tol && iter <= maxit

            try
                % Solve for sellers' entry threshold
                u0starnew       = v0starfunc(cL+eSo,eBo_posr_init,cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,screeningvalue,PAR);

                % letting lambdastar be defined as value that justifies u0star
                lambdastarposr = fminbnd(@(lambda) (u0starnew-v0starfunc(cL+eSo,eBo_posr_init,cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,screeningvalue,PAR,lambda)).^2,1,lambdastarposr+4);

                % Updating: 
                % seller parameters
                count = 1;
                err_count = 1;
                while (err_count == 1 || sum(thetas<ub)<length(thetas)) && count < 1001 
                    try
                        [thetasnew,fval]   = fmincon(@(theta) sellerLL(theta,PAR,data_posr,u0starnew,screeningvalue), startpts(count,:),[],[],[],[],PAR.lb,PAR.ub,[],PAR.opts);
                        err_count = 0;
                    catch 
                       count = count + 1;
                    end
                    count = count + 1;
                end

                % Update expected surplus matrices and resulting CDF
                [pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec]  = piBS_MC(cB,cS,thetab,thetasnew,PAR);
                CDFu0starnew = ggd_cdf(u0starnew,thetasnew);

                % Opportunity cost of time sellers
                eSonew      = eSofunc(screeningvalue,CDFu0starnew,lambdastarposr,cL,cR,thetasnew,pdfmat,pismat,v0vec,PAR);

                % To write to file
                output = [iter,CDFu0starnew,u0starnew,lambdastarposr,eSonew,thetasnew,fval];
                NPLiterations(iter+1,:)=output;
                
                 % Update diff 
                differences(1)=CDFu0starnew - CDFu0star;
                differences(2)=eSonew-eSo;
                differences(3:5)=thetasnew-thetas;
                diff = max(abs(differences))
                iter = iter + 1;
                CDFu0star = CDFu0starnew;
                u0star = u0starnew;
                eSo = eSonew;
                thetas = thetasnew;

            catch
                % Function values at interval endpoints must be finite and
                % real in v0* = didnt converge
                warning('No convergence. In v0starfunc: Stopped iterating because: Function values at interval endpoints must be finite and real.')
                iter = iter+1;
            end
        end
                
        % update eBo and eSo based on new u0star, original MLE estimate of lambdastarposr
        eBo_posr = eBofunc(screeningvalue,u0star,lambdastarposr_init,pibmat,pib0vec,pdfmat,v0vec,0,PAR); % second-last input: zeroreserve=0
        eSoMLE   = eSofunc(screeningvalue,CDFu0star,lambdastarposr_init,cL,cR,thetas,pdfmat,pismat,v0vec,PAR); 
        
     
        
        %%%% RESULTS %%%
    
        % Output parameters
        RobustInter = [thetas_init,eBo_posr_init,eSo_init,u0star_init,CDFu0star_init,screeningvalue,PAR.sharenobidders,iter]';
        RobustPars  = [thetab,thetas,eBo_posr,eBo_nor,eSoMLE,eSo,lambdastarnor,lambdastarposr_init,lambdastarposr,u0star,CDFu0star]';
        
        % Output mean / var V platform
        CDFb = @(x) ggd_cdf(x,thetab);
        Udraws         = arrayfun(@(i) fminbnd(@(x) (CDFb(x)-PAR.v0probs_fit(i))^2,-2,9),1:length(PAR.v0probs_fit))'; % Large sample of bidder conditional values
        meanVplatform = mean(Udraws);
        varVplatform  = var(Udraws);
        
        % Output mean / var V0 platform
        CDFs = @(x) ggd_cdf(x,thetas);
        U0draws         = arrayfun(@(i) fminbnd(@(x) ((CDFs(x)-CDFs(screeningvalue))./(CDFs(u0star)-CDFs(screeningvalue))-PAR.v0probs_fit(i))^2,screeningvalue,u0star),1:length(PAR.v0probs_fit))'; % Large sample of seller conditional values
        meanV0platform = mean(U0draws);
        varV0platform  = var(U0draws);
        
        RobustValues = [meanVplatform;varVplatform;meanV0platform;varV0platform];
        
        % Network effects adding one T / M in posr auctions
        lambdastarposr = lambdastarposr_init;
        T              = sum(data.hasreserve);
        M              = T*lambdastarposr ;
        changegrid     = [0,10];
    
        % Value of additional sellers / bidders (network effects) 
        T_grid             = changegrid+T;
        M_grid             = changegrid+M;
        v0bar_grid         = arrayfun(@(i) fsolve(@(x) (T/CDFu0star)*CDFs(x) - T_grid(i), u0star),1:length(T_grid));
        Pib_nsgrid         = arrayfun(@(i) Pibidder(eBo_posr,screeningvalue,v0bar_grid(i),M/T_grid(i),pibmat,pib0vec,pdfmat,v0vec,0,PAR), 1:length(T_grid));
        Pib_nsgrid_nov0bar = arrayfun(@(i) Pibidder(eBo_posr,screeningvalue,u0star,M/T_grid(i),pibmat,pib0vec,pdfmat,v0vec,0,PAR), 1:length(T_grid));
        Pis_nbgrid_marg    = arrayfun(@(i) Piseller(screeningvalue,CDFu0star,CDFu0star,eSoMLE+cL,eBo_posr,cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,PAR,M_grid(i)/T),1:length(M_grid));
        Pis_nbgrid_med     = arrayfun(@(i) Piseller(screeningvalue,CDFu0star,0.5,eSoMLE+cL,eBo_posr,cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,PAR,M_grid(i)/T),1:length(M_grid));
        
        RobustNE = [Pib_nsgrid(2)-Pib_nsgrid(1);Pib_nsgrid_nov0bar(2)-Pib_nsgrid_nov0bar(1);Pis_nbgrid_marg(2)-Pis_nbgrid_marg(1);Pis_nbgrid_med(2)-Pis_nbgrid_med(1)];
    
        % Lemons effect - figure 4
        quantilegrid = (0.1:0.1:0.8);
        PAR.highvalue=0;
         
        % Baseline surplus
        Pis_v0grid                  = arrayfun(@(i) Piseller(screeningvalue,CDFu0star,quantilegrid(i),eSo+cL,eBo_posr,cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,PAR,lambdastarposr),1:length(quantilegrid));

        % Increase cL (+1)
        Pis_v0grid_plus1cL_noentry  = arrayfun(@(i) Piseller(screeningvalue,CDFu0star,quantilegrid(i),eSo+cL+1,eBo_posr,cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,PAR,lambdastarposr),1:length(quantilegrid));
        v0star_plus1cL              = v0starfunc(cL+eSo+1,eBo_posr,cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,screeningvalue,PAR);
        Pis_v0grid_plus1cL          = arrayfun(@(i) Piseller(screeningvalue,CDFs(v0star_plus1cL),quantilegrid(i),eSo+cL+1,eBo_posr,cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,PAR),1:length(quantilegrid));

        % difference
        Lemondiff = Pis_v0grid_plus1cL-Pis_v0grid;

        %% Save results in Results structure
        % Fees
        Results.cS = cS;
        Results.cB = cB;
        Results.cL = cL;
        Results.cR = cR;
        Results.cP = 0;
        Results.eBo_posr = eBo_posr;
        Results.eBo_nor = eBo_nor;
        Results.eSo = eSo;

        % Estimated value parameters
        Results.mub            = thetab(1);
        Results.kb             = thetab(2);
        Results.mus            = thetas(1);
        Results.ks             = thetas(2);
        Results.shapeb=thetab(3);
        Results.shapes=thetas(3);
        Results.screeningvalue = screeningvalue;
        Results.sharetoolowr   = deleted_negative_impliedlnV0; %informative for fit

        % Entry 
        Results.u0star         = u0star;         %eqm, not structural
        Results.CDFu0star      = CDFu0star;      %eqm, not structural
        Results.lambdastarnor  = lambdastarnor;  %eqm, not structural
        Results.lambdastarposr = lambdastarposr; %eqm, not structural
        Results.sharenobidders = PAR.sharenobidders;   % =0 unless estimated 

        % First-stage fit
        Results.R2           = LM.Rsquared.Adjusted;

        % Some platform stats
        Results.T            = length(data.hasreserve);
        Results.NS           = round(Results.T/Results.CDFu0star);
        Results.M            = round(sum(data.hasreserve==0)*Results.lambdastarnor + sum(data.hasreserve==1)*Results.lambdastarposr);
        Results.p_r0         = 1-mean(data.hasreserve);
        Results.LMcoef   = LMcoef;

        % selected sample used in estimation - to use in modelfit etc
        if PAR.highvalue==0 && PAR.bootstrap==0
            writetable(data,'Estdata.csv');
            writetable(data_posr,'Estdata_posr.csv');
        elseif PAR.highvalue==1 && PAR.bootstrap==0
            writetable(data,'Estdata_high.csv');
            writetable(data_posr,'Estdata_posr_high.csv');        
        end


 

catch
        warning('Non-convergence / unstable (drop from bootstrap')
        Results.cS = 0;
        Results.cB = 0;
        Results.cL = 0;
        Results.cR = 0;
        Results.cP = 0;
        Results.eBo_posr = 0;
        Results.eBo_nor = 0;
        Results.eSo = 0;

        % Estimated value parameters
        Results.mub            = 0;
        Results.kb             = 0;
        Results.mus            = 999;
        Results.ks             = 999;
        Results.shapeb=0;
        Results.shapes=0;
        Results.screeningvalue = 0;
        Results.sharetoolowr   = 0; 
        Results.u0star         = 0;        
        Results.CDFu0star      = 0;  
        Results.lambdastarnor  = 0;  
        Results.lambdastarposr = 0;
        Results.sharenobidders = 0;  
        Results.R2           = 0;
        Results.T            = 0;
        Results.NS           = 0;
        Results.M            = 0;
        Results.p_r0         = 0;
        Results.LMcoef   = 0;

 end
    
end
